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I. INTRODUCTION 



We study the effects of two dynamical sea quarks on the spectrum of heavy quarkonia. Within the 
non-relativistic approach to Lattice QCD we found sizeable changes to the hyperfine splitting, but 
we could not observe any changes for the fine structure. We also investigated the scaling behaviour 
of our results for several different lattice spacings. 

o 
O 

00 . 

The experimental efforts to pin down the parameters of the Standard Model have been paralleled by intense 
, theoretical attempts to provide experimentalists with non-perturbative constraints from Quantum Chromodynamics 
(QCD). It is hoped that Lattice QCD will ultimately provide such important information. To this end it is crucial to 
£N| , understand and control the systematic errors in numerical calculations, which rely on extrapolations and interpolations 
to the physical point. This important task is particularly demanding for heavy quark phenomenology, where one has 
to describe accurately both the light and heavy quarks in the system. 

In particular the inclusion of light dynamical fermions in the gluon background is still a daunting task and requires 
the largest fraction of computational resources. In the past these restrictions led to the so-called quenched approx- 
imation, in which only valence quarks are allowed to propagate in a purely gluonic background, whereas the virtual 
creation of sea quarks is ignored. We have shown in a recent work JlJ that this results in systematic deviations in 
the lattice prediction of the light hadron spectrum from the observed experimental data. More recently it has also 
been found that the inclusion of two dynamical sea quarks has a significant effect on the light hadron spectrum and 
quark masses ^,||. This is of course analogous to QED, where the inclusion of all higher order effects, which could be 
made through perturbation theory, resulted in an impressive agreement with experimental observations. A distinctive 
aspect of QCD is that a proper non-perturbative treatment is in order so as to provide high-precision tests of this 
theory. Here we take this as our motivation to study heavy quarkonia in the presence (and absence) of dynamical sea 
quarks. 

Heavy quark systems have long been considered an ideal testing ground for QCD and they have triggered the 
development of static potential models [Q and heavy quark effective field theories 0. On the lattice, heavy quarks 
have frequently been studied using a non-relativistic approach to QCD (NRQCD fy) or a relativistic formulation 
promoted by the Fermilab group (t]]. Both provide effective descriptions to deal with large scale differences, which 
are difficult to accommodate on conventional lattices. NRQCD has been quite successful in reproducing the spin- 
independent spectrum of heavy quarkonia ^| owing to the fact that the quarks within such mesons move with small 
velocities v such that v 2 <C c 2 . 

As an effective field theory the predictive power of NRQCD relies on the control of higher dimensional operators, 
which have to be matched to the non-relativistic expansion of QCD. This has been the subject of many previous 
studies |9|-|ll| . As a result of these works it seemed plausible that bottomonium states could be accurately described 
in the NRQCD approach, whereas the spin structure of charmonium is very sensitive to the higher order relativistic 
corrections |l^]13| ]. 

Within the NRQCD framework sea quark effects on the heavy quarkonium spectra have previously been studied 
at a lattice spacing of a w 0.1 fm using the Kogut-Susskind |L4| or the Wilson |l5|] action for sea quarks. Here we 
present results for three lattice spacings in the range a ~ 0.2 — 0.1 fm, paying particular attention to the dependence 
on the sea quark mass and scaling properties. Our gauge configurations are generated with a renormalization-group 
(RG)-improved gluon action p6[ and a tadpole-improved clover quark action Jl7t for two dynamical flavours [Q. 
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Some measurements are also made on quenched configurations generated with the same gluon action for making 
direct comparisons of dynamical and quenched results. 

In Section |l| we introduce the actions which we use in our calculation. In Section III we give the details of our 
simulation, meson operators and fitting methods. Our results are discussed in Section IV and Section [y] concludes 
this paper. 



II. ACTIONS 



A. Glue: RG-Improved Action 



Since there is no unique discretisation of the continuum gluon action one can employ a set of operators to cancel 
some of the discretisation errors in the lattice action. The most common choice is to simply add a rectangular 1x2 
plaquette, Tri? M „, to the standard Wilson action, TrP^u, 

S s(9 2 ) = ^{coTriV+ciTWV} , (1) 

where Tr denotes the trace over all indices and cq + 8ci = 1. All such choices have the same continuum limit, but 
different discretisation errors. Here we adopt a prescription which is motivated by an RG-analysis of the pure gauge 
theory [c\ — —0.331 Jl^|). This has proven to be a suitable choice compared to, say, the Symanzik-improved action 
(ci = —1/12), for reducing scaling violation on coarse lattices. Instead of the coupling constant squared, g 2 , we often 
quote the value of [3 = 6/g 2 . 



B. Light Quarks: Clover Action 



The standard discretisation of the fermion action removes the doublers at the expense of 0(a) discretisation errors. 



It is possible to remove these errors by adding a single operator (the clover term) as first suggested in |17|: 



S q (g 2 ,m q ) = qQq = q{/^ + m q ) q + ar q A 2 q - c sw (g 2 ) ar -j-q a^F^ q . (2) 

Here the second term removes the doublers a la Wilson and the last is to reduce the resulting O(a) errors. We choose 
r = 1 and c sw = (1 — 0.1402 g 2 )" 3 ^ 4 . For the latter we follow the tadpole prescription of |L8]], which has been shown 
to improve the convergence of lattice perturbation theory significantly. Our choice is based upon the pcrturbative 
plaquette values as determined in |jqj . To one-loop order our choice differs from the correct value M] omv by 0.008g 2 . 
Hence we expect only small 0(aa) scaling violations due to radiative corrections from the clover action, and 0(a 2 ) 
errors from the gluon action. 

In our simulations we work with two flavours of degenerate quarks of a common mass: m q = m u = m c i- For further 
reference, it is customary to replace the bare quark mass by the hopping parameter: n = l/2(m q a + 4). Since the 
direct simulation of realistic light Wilson quarks is not feasible on present-day computers we study the spectrum at 
a sequence of different k. 



C. Heavy Quarks: NRQCD 



With the above actions we generated full QCD dynamical configurations on lattices of about 2.5 fm in spatial extent 
and lattice spacings ranging from approximately 0.1 to 0.2 fm. Such lattices are particularly suited to study light 
quark physics which is determined by a single scale: Aqcd ~ 200 MeV. However, for systems containing slow-moving 
and heavy quarks we have to adjust the theoretical description to take into account all the non-relativistic scales: 
mass (wq), momentum (rngw) and kinetic energy (rriQV 2 ). In this work we implement the NRQCD formulation 
to propagate the heavy quarks in a given gluon background. This approach has met with considerable success for 
b-quarks P,pT[. Also charm quarks have previously been studied in this framework JlJJl^]. 

Whereas the relativistic boundary value problem requires several iterations to determine the quark propagator, the 
NRQCD approach has the numerical advantage to solve the two-spinor theory as a much simpler initial value problem. 
The forward propagation of the source vector, S*(x), is described by: 
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-c ^ + c ° 2A(4) (4) 

The improved lattice operators Aj , E and B are defined as in f| . Other discretisations of NRQCD have been suggested 
in the past but they differ only at higher order in the lattice spacing. Here we follow [|ll|jl5| and employ a formulation 
which includes all spin terms up to 0(mv 6 ) in non-relativistic expansion of QCD. On the coarsest lattice we checked 
explicitly that our equation (|^) gives the same hyperfine splitting as from the asymmetric version employed in 
The parameter n was introduced to stabilise the evolution equation against high momentum modes. This is standard 
in such diffusive problems, but one should keep in mind that a change in n will have to be accompanied by a change 
in rriQ to simulate the same physical system. Alternatively one could decrease the temporal lattice spacing to prevent 
the high momentum modes from blowing up pfJ. 

For a single quark source at point y we have S(x) = i5^ 3 ^(x, y), but we also propagate extended objects with the 
same evolution equation (||). The operator Hq is the leading kinetic term and SH contains the relativistic corrections. 
The last two terms in 8H are present to correct for lattice spacing errors in temporal and spatial derivatives. For the 
derivatives we use the improved operators defined in || and we also replace the standard discretized gauge field 
by 

^ = §^-^(^(4^* + ^^ • (5) 

With this prescription we aimed to achieve an accuracy of 0(a 4 ) for the heavy quark sector. Of course we also 
expect terms of 0(aa 2 ) due to radiative corrections to this leading order result. In principle, we have to determine 
all coefficients in (^) by (perturbative or non-perturbative) matching to relativistic continuum QCD. Just as for the 
light quark sector we rely on a mean-field treatment of all gauge links to account for the leading radiative corrections. 
However, there is no unique prescription for such an improvement and several different schemes have been employed 
in the past. More recently it has been suggested that the average link variable in the Landau gauge should be less 
sensitive to radiative corrections since the gauge fields in Landau gauge have less UV fluctuations j2lj . Here we adopt 
this view and divide all links by the appropriate tadpole coefficient at each value of (/3, k): 

U M (x) — ► U^xj/uoL , u 0L = -{TrUft) lg ■ (6) 

An alternative and gauge-invariant implementation of mean-field improvement that has frequently been used in the 
past forces the average plaquette, P^ v , to unity: 

U^x) -> U^/uqp , u P = -(TrP^) 1 / 4 . (7) 

In some selected cases we have also implemented this method to estimate the effect of unknown radiative corrections 
to the renormalisation coefficients, Cj. In all applications of (|4|) we set them to their tree-level value 1. We denote as 
0(mv 6 , a 2 ) the evolution equation which includes all spin-dependent operators up to 0(mv e ) and where all operators 
have been improved to reduce the 0(a 2 ) errors. 
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III. SIMULATION 



A. Updates, Trajectories and Autocorrelations 

The gauge field configurations with two dynamical sea quarks used for the present study were generated on the 
CP-PACS supercomputer at the Center for Computational Physics, University of Tsukuba. They can be classified by 
two parameters, (/3, k), which determine the lattice spacing and the sea quark mass. A standard hybrid Monte Carlo 
algorithm is used to incorporate the effects of the fermion determinant. For the matrix inversion we implemented 
the BiCGStab algorithm. To reduce the autocorrelations between separate measurements we only used every fifth or 
tenth trajectory and binned the remaining data until the statistical error was independent of the bin size. In Tab. || 
we list the number of trajectories we generated for each set of couplings along with the other simulation parameters 
and the actual number of configurations we used in this study; for more details see [Q. The subsequent determination 
of the quarkonium spectrum is subject of this work. 

Since there has been no previous study of heavy quarkonia using the RG action for the gluon sector, we also 
supplemented our calculation by a comparative quenched calculation. The coupling constant [3 of these quenched 
configurations were chosen so that the string tension of the static quark- antiquark potential matches that of one of 
the dynamical runs. The parameters of these runs are also given in Tab. 1. 



B. Meson Operators 

To extract meson masses we calculate two-point functions of operators with appropriate quantum numbers. In a 
non-relativistic setting gauge-invariant meson operators can be constructed from the two-spinors x ( x )i ^(u) an d a 
Wilson line, W(x, y) : y (x)W(x, y)^>(y). The construction of meson states with definite J PC from those fundamental 
operators is standard [22] (on the lattice J labels the irreducible representations of the octahedral group (J = 
A 1 ,A 2 ,T 1 ,T 2 ,E)). 

Since here we are only interested in S and P-states it is sufficient to consider x (x) and x^ ( x ) (x) , respec- 
tively. The corresponding spin-triplets can be constructed by inserting the Pauli matrices into those bilinears. We 
also sum over different polarisations to increase the statistics. 

The overlap of those simplistic operators with the state of interest can be further improved upon. One way is to 
employ wavefunctions, which try to model the ground state more accurately. This requires assumptions about the 
underlying potential and gauge-fixed configurations. We decided to use a gauge-invariant smearing technique, which 
regulates the extent of the lattice operator, with a single parameter e: 

xHx) 6 -^x f (x) O (1-eA 2 ) 10 9(x) . (8) 

For computational ease wc limited this procedure to 10 smearing iterations and implemented it only at the source. 
From such operators we obtain the meson correlators as a Monte Carlo average over all configurations 

C e (x,y) = <tr [&(x,y) 6 G(x,y) (1-eA 2 ) 10 6* ]) , (9) 

where tr denotes contraction over all internal degrees of freedom. 

For the smeared propagator we solve (j^) with <S*(x, y) = <5(x, y) (1 — eA 2 ) 10 . We fix the origin at some (arbitrary) 
lattice point y = (y, 0). This creates a meson state with all possible momenta. In practice we employed up to 8 spatial 
sources on every configuration. This is permissible since heavy quarkonia are small compared to the lattice extent of 
about 2.5 fm. We explicitly checked the independence of such measurements by binning. At the sink, x = (x, t), we 
perform a Fourier transformation to project the correlator onto a given momentum eigenstate: 

C e , t (p) = ^C e (x,t)exp(-*px) . (10) 

X 

In the trivial case of zero momentum this amounts to simply summing over all spatial x. 
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C. Fitting 



Since there is no backward propagation of the heavy quark in our framework, we can fit the meson propagators to 
the exponential form: 

VttfaEk) = J2 ai(a, p, e) e ~ E ^ * . (11) 
»=i 

This is the theoretical prediction for a multi-exponential decay of a state with momentum p and quantum numbers 
a along Euclidean time t. Different choices for the smearing parameter e will result in different overlaps with the 
ground state and the higher excited states of the same quantum number. In practice it is difficult to project directly 
onto a given state, so we chose to extract the ground state from multi-exponential correlated fits. In some cases we 
were also able to extract the first excited states reliably from our data. The simplest way to visualise our data is by 
means of effective mass plots, which are expected to display a plateau for long Euclidean times. In Fig. [I] we show a 
representative plot for one set of simulation parameters. 

Correlations between different times, t, and different smearing radii, e, of the meson propagator C £j t are taken into 
account by employing the full covariance matrix for the ^-minimisation: 

X 2 (tt) ee (Or - ifrM) r- 1 (C s - y s (n)), ■ (12) 

r,s — 1 

Here, in order to ease the notation, we introduced multi- indices (r = [e, <]) for the data points and 7r = ja-j, Ek] for the 
parameters. 

Our statistical ensembles are large enough to determine the covariance matrix, r~ s , with sufficient accuracy. 
Therefore the inversion of this matrix did not present a problem. For the spin splittings we applied ( |ll| ) also to 
the ratio of two propagators, C (a.\) / C (a.2) ■ In this way we utilized correlations between states of different quantum 
numbers to obtain improved estimates for their energy difference. 

Statistical fluctuations in the data cause fluctuations in the fit results determined by (|l2|). We estimate the 
covariance matrix, Aj~i, of the fitted parameters, 7Tfc, from the inverse of (d 2 x 2 ) / '{dukd^i). We have checked these 
errors against bootstrap errors which gave consistent results. We also require consistent results as we change the 
fit ranges (t m i n ,t max ) or the number of exponentials to be fitted. This is illustrated in Fig . H , where we show the 
twin-plot for the S-state of Fig. [jj The goodness of each fit is quantified by its Q-value p3fl and we demand an 
acceptable fit to have Q > 0.1. Finally we subjected those results to a binning analysis, which takes into account 
autocorrelations of the same measurement at different times in the update process. 



IV. RESULTS 



We now present our new results for elements of the spectrum of heavy quarkonia. Our data from two quenched 
lattices is given in Tab. || and the dynamical data is collected in Tabs. Ill to VI. For notational ease we use 
dimensionlcss lattice quantities throughout the remainder of this paper, unless stated otherwise. To convert the 
lattice predictions into dimensionful quantities we take the experimental IP — IS splitting to set the scale. 



A. Heavy Mass Dependence and Kinetic Mass 

For a given value of the gauge coupling (lattice spacing) and the mass of the two degenerate sea quarks there is 
only one remaining parameter to choose - the mass of the heavy valence quark. On the lattice we are free to simulate 
every arbitrary value, but in order to obtain physical results we tune the ratio M^i n /(1P — IS) of the kinetic mass 
of a quarkonium state and the IP — IS mass splitting to its experimental value. The determination of 15* and IP 



masses has already been described in Sec. Ill C . The kinetic mass, M^ in , is defined through the dispersion relation of 
the quarkonium state: 

E(p) - E(0) = + . . . , p = ^(ni,n a ,n 3 ) . (13) 
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For each heavy quark mass, tuq, we project the 1 So-state and the 3 Si-state onto 5 different momentum eigenstates: 
(ni,n 2 ,n 3 ) = (0, 0, 0); (f , 0, 0); (f , 1 , 0); (f , 1, 1 ) and (2,0,0). We obtain E(p) - E(0) from ratio fits and determine 
the kinetic mass by fitting the data to the dispersion relation. To this end we have also included higher terms in the 
expansion of ([l|) and find consistent results for Afki n . However, such fits tend to have larger errors and the coefficient 
of p 4 is not well determined. For better accuracy we normally restrict ourselves to a linear fit in p 2 for the lowest 
four momenta. An example of this procedure is given in Fig. ||. We have plotted the fitted values of Mkin against 
the heavy quark mass in FigM Large discretisation errors can be seen for almost all masses, but once we include all 
0(a 2 ) correction terms in (0)^ the discrepancy due to higher order relativistic corrections is small. Comparing the 
relative changes in Fig. || due to 0(a 2 ) improvement at different momentum scales, we can also estimate the size of 
higher order corrections and we expect them to be small. 

In this paper we tuned the bare quark mass on all our lattices n), so as to reproduce the experimental value of 
the mass of T (Mki n = 9.46 GeV). In some selected cases we interpolated the spectrum to this physical point. 



B. Scale Determination and IP IS Splitting 

It has been noticed in the past that the tuning of tuq can be done efficiently since the spin-averaged splitting is 
not very sensitive to the quark mass parameter. However, with our newly achieved accuracy we could also resolve a 
slight mass dependence of IP — IS in the range from charmonium to the bottomonium system as shown in Fig. |^. 
The experimental values for the IP — 15 splitting show a 4% decrease when going from charmonium (458 MeV) to 
bottomonium (440 MeV), which should be compared to a reduction of about 10% from our simulation at Nf = 2 and 
an unphysical sea quark mass. This larger change is related to the running of the strong coupling between the two 
scales, which still does not fully match the running coupling in nature. It is expected that the modified short-range 
potential will result in a ratio (IP — lS) ce /(lP — lS)^ bigger than in experiment. 

While the heavy quark mass can be tuned to its physical value as described in the previous section, this is not 
possible for the light quark mass and one has to rely on extrapolations to realistic quark masses, where the ratio 
m^/nip equals the experimental value. Here we are mainly interested in the behaviour of physical quantities as we 
approach the chiral limit. We take m 2 as a measure of the light quark mass and extrapolate quadratically in this 
parameter. This is a common procedure but we will demonstrate below that the physical dependence on the sea 
quark mass may indeed be difficult to disentangle from unphysical scaling violations. In taking the naive chiral limit 
we hope to account for at least a fraction of the spectral changes towards smaller sea quark masses. At /3 = 2.10 we 
only have results from two values of k and take a linear estimate for the chiral limit. The chiral behaviour of the 
IP — IS splitting is shown in Fig. || for all values of f3 in our study. 

In quenched simulations there exist uncertainties when setting the scale from different physical quantities. We 
expect these uncertainties to be reduced in our simulations incorporating two light dynamical flavours. To examine 
this point we compare in Fig. (?] our results for IP — IS splitting with the data for m p as a representative example 
from the light quark sector. If it were not for quenching effects and lattice spacing artefacts, one would expect the 
ratio m p /(lP — IS") to equal its experimental value. 

It is encouraging to see that the dynamical calculations are always and significantly closer to the experimental value 
of 1.75 than the corresponding quenched simulations. This demonstrates the importance of dynamical over quenched 
simulations. The scaling violations in this ratio do not fully cancel, however; we observe a 10% shift in m p /(lP — 15) 
over a « 0.2 — 0.1 fm. Keeping in mind that we are working on rather coarse lattices with a>0.1 fm, the remaining 
scaling violations are perhaps not too surprising. 

Looking at the ratio (IP — 15)/Aqcd it is clear that our data does not satisfy the strict criterion of asymptotic 
scaling, see Fig. ||. In this plot we determined Aqcd from the 2-loop formula in the MS scheme, 

. fab \ { - bl/2 ^ ( 2tt\/ bl-b 2 b Q \ 

Aqcd ^U^ ex H~w v 1+a_ ^rJ ' (14) 

where the MS coupling constant a = aj^g(ir / a) is estimated with a tadpole-corrected one- loop relation defined by 

1 _ (3.648P - 2.648P) 



a JTs( 7T / a ) a o 



4^(0.0589 + 0.0218^/) . (15) 



Here a = ,9 2 /47r denotes the bare coupling, and P and R are, respectively, lxl and 1x2 Wilson loops normalized 
to unity for U^x) = 1. 
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Within the effective approach of NRQCD, we cannot extrapolate such scaling violations away and it is crucial to 
find other ratios in which the scaling violations cancel each other already at finite lattice spacing. In Fig. ^ we show 
a test of this nature for the string tension, which shows a better scaling. Here we plot as open symbols the data 
obtained from 4 different sea quark masses. At (3 = 2.20 (2.10) we only measured the lightest (and heaviest) sea-quark 
mass, corresponding to m T /m p « 0.60 (0.80). This figure also suggests that the string tension, in units of the IP — IS 1 
mass splitting, is smaller for 2-flavour QCD when compared to the quenched (Nf = 0) theory. 



C. Hyperfine Splitting 

Quenching effects are also expected to show up in short-range quantities, since they are particularly sensitive to 
the shape of the QCD potential. In H[m| this difference has been demonstrated explicitly by observing a change in 
the Coulomb coefficient of the static potential. In the context of heavy quarkonia, the hyperfine splitting is such a 
UV-sensitive quantity which should be particularly susceptible to changes in the number of flavours and the sea quark 
mass. The prediction from potential models is 



In our study this is the most accurately measured quantity and it is clearly very sensitive to the value of the heavy 
quark mass, see Fig. [l^. As has been noticed previously, higher order relativistic and radiative corrections are 
equally important for precision measurements of the hyperfine splitting in bottomonium |ll|,|l5| and even more so for 
charmonium fl3|| . Here we employ 0(mv e , a 2 ) as the standard accuracy throughout this paper. 

The equation (|l^) involves a direct dependence on both the strong coupling and the wavefunction at the origin, 
which makes the hyperfine splitting an ideal quantity to study quenching effects. Here we also observe a clear rise of 
the hyperfine splitting as we decrease the sea quark mass, see Fig. [ll]. 

In Fig. |l2] we collected all our dynamical results for the hyperfine splitting over the range of 0.1 — 0.2 fm. Here we 
plotted the data from each sea-quark mass as open symbols and used the experimental IP — IS* splitting to convert 
lattice data into MeV. One should keep in mind that these points correspond to unphysical bottomonium in a world 
of different sea quark masses. We also plot as filled symbols the results of our naive chiral extrapolation as described 
in the previous section. 

At around 0.10 fm we notice a very good agreement with the only previous calculation jl^]. These authors have 
performed a dynamical simulation at a single lattice spacing using Wilson glue and unimproved sea-quarks. For the 
bottom quarks they used an NRQCD formulation with the same accuracy, 0(mv 6 , a 2 ), as in this study. 

An unpleasant feature with our results in Fig. |l^ is lack of scaling; for both the full and quenched case we find 
scaling violations of about 100 MeV/fm for the hyperfine splitting. Nonetheless, we do find several indications for sea 
quark effects in our results. First we notice that, if it were not for sea quark effects, then all points in Fig. [l~2| would 
lie on a universal curve which is not the case. This is a strong indication that for this quantity we have to expect 
effects of the order of 3-5 MeV when going from zero to two flavour QCD. 

To substantiate this observation we make a direct comparison of quenched and dynamical calculations at the same 
lattice spacing of 0.14 fm in Fig. [l^. For the 3 Si — 1 Sq splitting replotted from Fig. |l2|, a clear increase of around 5 
MeV (20%) represents more than a 5cr-effect, which reflects the accuracy in our determination of this quantity. On 
the other hand, the hyperfine splitting in P-states is reduced as we approach a more realistic description of QCD. 
Within potential models, states with L ^ are not sensitive to the contact term of the spin-spin interaction. However 
the perturbative expression for a higher order radiative corrections |p6| gives a 3 P — 1 Pi splitting opposite in sign to 
our values. Experimentally, the spin-triplet states are well established, but 1 S'o and 1 P\ have yet to be confirmed for 
bottomonium. 



We comment that our data for charmonium (Tab. VI) also indicates a rise in the hyperfine splitting towards the 
chiral limit. It is, however, apparent that such a rise can not explain the discrepancy between the NRQCD predictions 
and the experimentally observed spin structure. We confirm an earlier observation |l3| ] that the velocity expansion is 

0.3. 



D. Fine Structure 



In the continuum, the fine structure in quarkonia is due to the different ways in which the spin can couple to the 
orbital angular momentum of the bound state. In our approach, the spin-orbit term and the tensor term of potential 
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models can be traced back to the C3 term in (Q). A correct description of the fine structure will therefore require a 
proper determination of C3 and the corrections to this term. 

On the lattice we have also additional splittings with no continuum analogue. For example, the 3 P2E — z Pn 
splitting is known to be a pure discretisation error since the lattice breaks the rotational invariance of the continuum 
and causes the J = 2 tensor to split into two irreducible representations of the orthogonal group: T2 and E. Indeed, 
for both the dynamical calculations as well as the quenched data, we observe a significant reduction of this splitting 
when the lattice spacing is decreased; the splitting diminishes from about 1.5 MeV at a « 0.2 fm to 0.5 MeV at 
a ~ 0.1 fm for the dynamical case. 

In Fig. [[4] we show our results for the fine structure. For and 3 Pi we observe no clear dependence on the sea 
quark mass. This is not totally unexpected since P-state wavefunctions vanish at the origin and should not be as 
strongly dependent on changes in the UV-physics. In any case such small changes would be difficult to resolve within 
our statistical errors. 

From Fig. |IJ we can also see a better scaling behaviour of the P states, apart perhaps from the 3 Po, where scaling 
violations still obscure the chiral behaviour. The latter has J = and therefore we may expect that for this state 
restoration of rotational invariance is particularly important. 

On our finer lattices we observe an increase of the 3 Po — 3 P splitting, closer towards the experimental value of 
—40 MeV. We take this as an indication that a better control of the lattice spacing errors and radiative corrections 
is necessary to reproduce this quantity in NRQCD lattice calculations. The other splittings, 3 P2 — 3 P and 3 Pi — 3 P, 
deviate only by a few MeV from their experimental values of 12 MeV and —8 MeV, which could be due to missing 
dynamical flavours (Nf — 3), higher order relativistic effects and radiative corrections. 

We take the fine structure ratio, Rf s = ( 3 P2 — 3 Pi)/( 3 Pi — 3 Pq), as a particularly sensitive quantity to measure 
the internal consistency of the P-triplet structure. This quantity should be less sensitive to radiative corrections of 
the NRQCD coefficients away from their tree-level values. Previous NRQCD calculations had measured this quantity 
to be much larger than 1, compared to the experimental value of 0.66(4). We believe that this discrepancy is due to 
lattice spacing artefacts as it is very sensitive to the implementation of 0{a 2 ) improvement in the NRQCD formalism. 
It is encouraging to see that this value is further reduced on our finer lattices, see Fig. [[J. Notably, we do not observe 
any difference between our dynamical results and the quenched data. 

E. 2S - IS Splitting 

Another spectroscopic quantity which has attracted much attention is the 2S — IS splitting, since it should also 
be sensitive to the short-range potential. On conventional lattices such higher excitations are difficult to resolve and 
requires delicate tuning to minimise the mixing of the 2S with the ground state. Given our rather coarse lattices we 
did not attempt to perform a systematic study of this quantity, but in the context of this section it is important to 
notice that we do not observe any chiral dependence of the ratio, R2S = (2S 1 — IS) /(IP — IS). 

In Fig. [R3| we compiled representative data from other groups |l5| , p5|] along with our new results from the RG-action. 
Within the large errors we cannot resolve a discrepancy between the experiment and the lattice data. This result is 
in contrast to the previous determinations of this quantity which claim to see deviations due to missing sea quarks 

EM- 

Observing such deviations is certainly plausible as this ratio is thought to be sensitive to differences in the underlying 
short-range potential. However, for the same reason we should also expect large scaling violations. Interestingly, on our 
coarsest dynamical lattices we even observe smaller values of i?2S, which we take as an indication of large discretisation 
errors. Apart from this very coarse lattice data, we cannot resolve either scaling violations or quenching effects. We 
feel that it requires a much better resolution of the higher excited states, which is hard to achieve on isotropic lattices. 
Future lattice studies will need optimised meson operators or finer temporal discretisations to observe these effects. 

V. CONCLUSION 

We have demonstrated that dynamical sea quarks have a significant effect on the spectrum of heavy quarkonia. 
Namely the hyperfine splitting, 3 Si — So, is raised by almost 20% when going from zero to two dynamical flavours. 
The efficiency of the NRQCD approach has played an important role to establish such effects, but the numerical 
simplicity of this approach is offset by additional systematic errors, which have to be controlled. The sensitivity of the 
spin-structure to relativistic, 0(mv G ), and radiative, 0(a), corrections was well known before we started this work. 
Here we demonstrated that quenching errors are equally important for precision measurements of the spectrum of 
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heavy quarkonia. Perhaps more worrying are scaling violations, which we could resolve in many quantities. Without 
a proper control of lattice spacing artefacts it is not possible to make predictions for such UV-sensitive quantities as 
the hyperfine splitting on the lattices we used here. 

While the lattice predictions for 3 Pi — 3 P and 3 P2 — 3 P agree well with their experimental values, the determination 
of the 3 Pq is more problematic and we still observed large deviations from the experimental value when the IP — 15 
splitting is used to set the scale. Clearly much work remains to be done to reduce the systematic errors in heavy 
quark physics to the same degree as the statistical ones. We feel that this may be difficult to achieve within the 
non-relativistic framework. A better description of the NRQCD coefficients or a relativistic treatment is in order to 
describe the spin structure in quarkonia accurately. From this work it is apparent that full QCD simulations are also 
necessary to achieve such a goal. 

For less accurate quantities such as 25 — 15 it is more difficult to reach a similar conclusion and we leave a decisive 
observation of both scaling violations and sea quark effects to future studies with refined methods. 
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0.8805090(44) 


0.850170(25) 
0.854604(20) 


2000/192 
2000/192 


2.20, (24 3 x 48) 


1.44 


0.1368 


1.95, 2.21 


0.8878887(29) 


0.866357(23) 


2000/128 








quenched simulation, Nf 


= 






(3, (L 3 x T) 








UOL 


updates/cfgs. 


2.187, (16 3 
2.281, (16 3 


x 32) 
x 32) 




3.70, 4.00 
3.20, 3.40 


0.8772362(22) 
0.8855537(18) 


0.831789(55) 
0.847829(41) 


20000/200 
20000/200 



TABLE I. Simulation parameters for this study. The last column gives the number of configurations employed for this work. 



The quenched runs are made at /3 = 2.187 and 2.281 so that the string tension matches with those of the Nf — 2 runs at 
{13, k) = (1.95,0.1390) and at (1.95, k c ) on a 16 3 x 32 lattice. 



P 


2.187 


2.187 


2.281 


2.281 


M b 


3.70 


4.00 


3.20 


3.40 


M kin [GeV] 


9.04(23) 


9.95(27) 


9.110(94) 


9.77(10) 


a(lP - IS) [fm] 


0.1653(23) 


0.1637(15) 


0.1423(12) 


0.1400(12) 


R2S 


1.50(25) 


1.65(59) 


1.299(98) 


1.26(11) 


i S 1 - 'So [MeV] 


23.12(34) 


21.68(22) 


24.88(25) 


23.91(26) 


*P - 1 P 1 [MeV] 


4.03(61) 


3.79(61) 


4.97(35) 


4.87(36) 


3 P - S P [MeV] 


20.4(1.3) 


19.5(1.4) 


29.8(1.0) 


28.1(1.2) 


3 P - 3 Pi [MeV] 


6.73(65) 


6.46(74) 


9.07(58) 


8.93(59) 


3 P 2 - 3 P [MeV] 


6.83(56) 


6.44(65) 


9.26(54) 


9.23(51) 


3 P 2 e - 3 P 2 t [MeV] 


1.67(45) 


1.43(28) 


0.91(36) 


0.87(35) 


a P 2 - 3 Pi [MeV] 


13.6(1.2) 


12.9(1.4) 


18.8(1.1) 


18.2(1.1) 


3 Pi - 3 P [MeV] 


13.22(76) 


12.73(82) 


20.53(62) 


19.39(65) 


Rfs 


1.03(11) 


1.01(13) 


0.917(58) 


0.937(64) 



TABLE II. Bottomonium spectrum from quenched run at f3 = 2.187 and (3 — 2.281. These results should be compared 
directly to the Nf = 2 data at the similar lattice spacing: ((3, k) = (1.95, 0.1390) and (1.95, k c ), respectively. We also illustrate 
the effects of little changes in the quark mass on the spectrum. The difference for the hyperfine splitting, 3 Si — 1 So, is most 
noticeable. For the other splittings we see indications of such a suppression as the mass is increased, but it is much less 
pronounced within the errors. Scaling violations can be observed in P2E — Ptt, as discussed in the main text. All the other 
splittings are suppressed on the coarser lattice. On the finer lattice, the ratio Rf s = ( 3 P 2 — 3 Pi)/( 3 Pi — 3 Po) is closer to its 
experimental value: 0.66(4). 
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K 


0.1409 


0.1430 


0.1445 


0.1464 






0.80599(75) 


0.7531(13) 


0.6959(20) 


0.5480(45) 




M b 


5.87(5) 


5.85 


5.61 


5.10(5) 




M kin [GeV] 


9.46(12) 


9.43(10) 


9.530(80) 


9.46(11) 




a(lP — IS) [fm] 


0.2787(25) 


0.2765(26) 


0.2611(19) 


0.2306(21) 


0.1987(32) 


"■(ppq) [f m ] 


0.2622(11) 


0.2560(16) 


0.2462(13) 


0.2246(18) 


0.2040(40) 


R 2 s 


1.157(25) 






1.143(52) 




3 Si - 'So [MeV] 


20.80(33) 


19.75(29) 


21.11(24) 


23.50(35) 


26.81(76) 


3 P - x Pi [MeV] 


3.75(24) 


3.91(18) 


3.64(54) 


3.40(56) 


2.82(1.03) 


3 P - 3 P [MeV] 


11.80(59) 


12.26(65) 


13.09(33) 


11.8(1.1) 


13.61(75) 


3 P - 3 Pi [MeV] 


5.71(33) 


5.20(30) 


5.78(15) 


4.57(64) 


5.64(12) 


3 P 2 - 3 P [MeV] 


6.08(31) 


5.97(33) 


6.23(16) 


5.51(64) 


6.12(39) 


3 P 2 b - 3 P 2 t [MeV] 


1.84(23) 


1.56(26) 


1.67(13) 


1.58(47) 


1.47(31) 


3 P 2 - 3 Pi [MeV] 


11.81(63) 


11.18(60) 


12.02(30) 


10.1(1.3) 


11.69(77) 


3 Pi - 3 P [MeV] 


6.02(31) 


7.19(42) 


7.33(23) 


6.91(63) 


8.13(46) 


Rfs 


1.96(14) 


1.55(12) 


1.640(66) 


1.46(23) 


1.44(13) 



TABLE III. Bottomonium results from (3 = 1.80. An error in the quark mass parameter indicates that we have interpolated 
to this value in order to reproduce the Bottomonium at any given (/?, k). Where this error is not given we accepted the value 
of the tuned quark mass. For the hyperfine splitting we could fit the data to a linear-plus-quadratic dependence in the quark 
mass, but for the fine structure the quadratic part was ill-determined and we resorted to linear or constant fits if their Q-value 
was acceptable. 
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0.1410 




m^/mp 


0.80484(89) 


0.7514(14) 


0.6884(15) 


0.5862(33) 




M b 


4.00 


3.80 


3.70 


3.40 




M kin [GeV] 


9.40(16) 


9.43(22) 


9.43(10) 


9.49(17) 




o(lP - IS) [fm] 


0.1767(13) 


0.1662(35) 


0.1586(15) 


0.1473(17) 


0.1341(48) 


a(ppg) [fm] 


0.1974(11) 


0.1860(12) 


0.1791(10) 


0.1625(13) 


0.1451(33) 


R 2S 




1.242(72) 


1.46(17) 


1.47(31) 




3 Si - 'So [MeV] 


25.11(49) 


25.72(81) 


26.07(38) 


27.85(60) 


30.07(1.58) 


3 P - x Pi [MeV] 


2.41(66) 


2.26(63) 


2.50(64) 


2.67(32) 


2.55(24) 


3 P - 3 P [MeV] 


18.4(1.8) 


20.0(1.8) 


21.5(1.7) 


23.1(1.7) 


24.9(5.5) 


3 P - 3 Pi [MeV] 


6.08(99) 


5.97(93) 


6.38(82) 


5.83(86) 


5.5(2.2) 


3 P 2 - 3 P [MeV] 


8.7(1.0) 


9.01(93) 


8.98(85) 


9.10(89) 


9.0(2.3) 


3 P 2 e - 3 P 2 t [MeV] 


2.05(15) 


1.50(13) 


1.41(20) 


1.22(74) 


1.33(29) 


3 P 2 - 3 Pi [MeV] 


14.8(2.0) 


15.1(1.8) 


15.5(1.6) 


14.8(1.7) 


14.2(4.4) 


3 Pi - 3 P [MeV] 


12.4(1.1) 


13.2(1.0) 


14.9(1.1) 


17.4(1.0) 


20.9(2.5) 




1.19(19) 


1.14(16) 


1.04(13) 


0.85(11) 


0.68(23) 



TABLE IV. Bottomonium results from j3 = 1.95. 
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TABLE V. 


Bottomonium results from f3 = 


2.10 and 2.20. 





CM) 


(1.80,0.1409) 


(1.80,0.1430) 


(1.80,0.1445) 


(1.80,0.1464) 


(1.95,0.1375) 


exp. 


iiIk /m p 


0.80599(75) 


0.7531(13) 


0.6959(20) 


0.5480(45) 


0.80484(89) 


0.18 


M b 


2.20 


2.10 


2.06 


1.77 


1.30(5) 




M kin [GeV] 


3.019(87) 


3.323(34) 


3.589(46) 


3.401(85) 


3.01(12) 


3.09688(4) 


a(lP - IS) [fm] 


0.2874(11) 


0.2758(14) 


0.2571(26) 


0.2388(53) 


0.1983(43) 




a(ppo) [fm] 


0.2622(11) 


0.2560(16) 


0.2462(13) 


0.2246(18) 


0.1974(11) 




i?2S 


1.378(60) 


1.29(10) 


1.557(95) 


2.02(34) 




1.3009(31) 


3 Si - [MeV] 


49.60(35) 


53.17(38) 


54.02(67) 


56.04(70) 


55.5(2.8) 


117(2) 


3 P - x Pi [MeV] 


3.66(23) 


2.86(86) 


3.25(41) 


1.71(55) 


4.0(2.0) 


-0.86(25) 


3 P - 3 P [MeV] 


26.48(54) 


31.52(74) 


26.6(3.2) 


31.4(2.1) 


38.8(4.1) 


110.2(1.0) 


3 P- 3 Pi [MeV] 


5.14(32) 


7.17(47) 


2.6(1.7) 


4.27(75) 


0.90(1.0) 


14.75(18) 


3 P 2 - 3 P [MeV] 


8.86(79) 


10.85(49) 


7.9(1.6) 


10.54(64) 


7.2(1.0) 


30.89(18) 


3 P 2 b - 3 P 2 t [MeV] 


2.45(18) 


3.00(46) 


2.06(30) 


2.18(66) 


1.50(78) 




3 P 2 - 3 Pi [MeV] 


13.12(49) 


18.15(76) 


10.4(3.2) 


14.8(1.3) 


9.3(4.0) 


45.64(18) 


3 Pi - 3 Po [MeV] 


21.17(33) 


24.40(49) 


24.23(99) 


29.25(84) 


34.8(5.0) 


95.4(1.0) 


Rfs 


0.620(25) 


0.744(34) 


0.43(13) 


0.507(47) 


0.27(12) 


0.4783(54) 



TABLE VI. Charmonium results. 
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FIG. 1. A representative effective mass plot for S- and P-states at (/?, k, 7tiq) — (1.80,0.1409,6.00). One can clearly observe 
a plateau for long enough times. For the S-states we employed 3 different smearings, which result in different overlaps with the 
ground state. 
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FIG. 2. Here we show consistent fit results for the S-state of Fig. [l]when plotted against the start of the fit range, train- We 
fixed tmax = 24 throughout. Different symbols denote different values for nat in the multi-exponential fit of Eq. O. 
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FIG. 3. These figures illustrate the tuning of the quark mass as described in the main text. On the left-hand side we show 
the i min -plots for the ratio fits of different momentum states with respect to the 3 Si at rest. We can perform two consistent 
fits up to p 2 (dashed line) and up to p 4 (solid line) in the dispersion relation, Eq. [l^. 
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FIG. 4. Quark mass dependence of Mkin- Here we can see sizeable discretisation errors for almost the whole range of quark 
masses between charm and bottom at (/3, «) = (1.80,0.1409). The implementation of 0(a 2 ) improvement in the NRQCD 
approach is clearly important on our lattices. In contrast, the sensitivity of A/kin to the relativistic correction terms is much 
smaller. 
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FIG. 5. Heavy quark mass dependence of the IP — IS splitting. We plot the (IP-IS) splitting against the heavy quark mass 
at (f3, k) = (1.95, 0.1375) and with two value of the stability parameter, n = 1,2. The vertical lines denote the regions of the 
charmonium and bottomonium system. 
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FIG. 6. Light quark mass dependence of IP — IS splitting in Bottomonium. We use quadratic fits in m% to extrapolate our 
data from four different sea quark masses to the chiral limit. For the two sea quark masses at (3 = 2.10 we show an estimate 
of the chiral limit by drawing straight lines. The single point at (/3, k) = (2.20, 0.1368) is shown for completeness. 
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FIG. 7. Here we show the ratio m p /(lP — IS), where scaling violations can be seen. In each case we use open symbols 
to denote data from dynamical calculations with different sea quark mass and full symbols to mark the chirally extrapolated 
values. Representative quenched results are also shown as full symbols. 



1G 



3.00 



2.80 



8 2.60 
w 

°- 2.40 



2.20 



2.00 
0.00 



0.04 



0.08 0.12 
a iP-is [ fm ] 



0.16 0.20 0.24 



FIG. 8. Asymptotic scaling. In this plot we take the chirally extrapolated values for IP — 15* and compare their scaling 
behaviour with respect to Aqcd- The latter is taken from 2- loop perturbation theory. 
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FIG. 9. In contrast to Figs. |^ and ^, we observe a better scaling for the ratio *J~o / (IP — IS) on our finer lattices. 
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FIG. 10. Mass dependence of hyperfine splitting. Here we show the strong mass-dependence of the hyperfine splitting, 
3 Si — 1 So, plotted against the inverse kinetic mass at (/3,k) = (1.95,0.1375). The vertical lines denote the regions of the 
bottomonium and charmonium system. This splitting is clearly very sensitive to the parameters of NRQCD. All data points 
are from updates with 0(mv 6 ,a 2 ) and n = 1,2 denotes different values of the stability parameter. 
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FIG. 11. Hyperfine splitting vs. m^. Here we collect the data for the the 3 Si — 1 S'o splitting in bottomonium from all values 
of (/3, k). A clear dependence on the sea quark mass can be seen. The linear-plus-quadratic fit curves are shown as solid lines. 
Here we used the IP — IS splitting to determine the lattice spacing. 
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FIG. 12. Scaling violations for hyperfine splitting. Open symbols correspond to runs with different sea quark mass. Filled 
symbols denote the dynamical data after chiral extrapolation and results with Nf = 0. We used IP — 15 splitting to determine 
the lattice spacing. 
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FIG. 13. Direct comparison of the bottomonium spin structure for quenched and full QCD at the same lattice spacing of 
a =i 0.14 fm. The Nf = 2 data is taken from the chiral limit of our measurements at /? = 1.95. 
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FIG. 14. Fine structure in Bottomonium. Here we plot (top to bottom) 3 Pi, 3 Pi and 3 Po relative to the spin averaged triplet 
state: 3 P = 1/9(5 3 P 2 + 3 3 Pi + 1 3 P )- The corresponding experimental values are: 13 MeV, -8 MeV and -40 MeV. 
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FIG. 15. Fine structure ratio in Bottomonium. Here the lattice data should be compared to the experimental value of 
0.66(4). It is apparent that there are still large underlying scaling violations, but no clear sea quark dependence. 
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FIG. 16. Here we show a scaling plot of the ratio R2S = (2S 1 — 15')/(1P — IS). It is apparent that one needs much smaller 
statistical errors to resolve any systematic effects. We show our results from different sea quark masses along with representative 
results from other collaborations. 
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